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LOADING - PLANE PROBLEM, PART I 



1. Introduction 



In Ref. 1 possible problem formulations for predicting 
the effects of cavitation on underwater shock propagation were 
examined. It was concluded that using either the displacement 
potential (a scalar) or the particle displacement (a vector) 
as the dependent variable could produce useful results. 

Ref. 2 gave a summary account of an unsuccessful effort to 
use ADINA for applications with axi symmetric geometry. 

The present report is a summary of progress to date in 
developing a capability for modeling a two-dimensional fluid 
region in which a structure is submerged. Examination of 
available codes disclosed no candidate allowing implementation 
of the displacement potential formulation, \vith cavitation, 
for the fluid and simultaneously permitting an appropriately 
coupled structural model. Accordingly, an ad hoc code for 
the fluid has been developed. A compatible structural code 
is being constructed. All codes referred to in this report 
are written in FORTRAN IV. 



2.1 Mesh Generator (MGNEW) 

For a plane problem with a cvlindric structure whose 
length is normal to the plane, the fluid region is the portion 
of the plane exterior to a circle (radius A) . Taking advantage 
of symmetry and truncating gives the rectangular region of 



2. Program Development 
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Fig. 1. The mesh is generated from a set of super - e lements 
as shown in Fig. 2. The transitional super-elements (those 
adjacent to the structure) each have two curved sides. The 
remaining elements are square. Note that, if each quadrant 
of the structural "hole" is treated as an element, the region 
is topologically equivalent to MU layers, each containing 
ML + MR squares. In the region of Fig. 2 the entry face is 
at x = MR*(0.9*A) and the exit face at x = ML*(0.9*A). The 
top face is at y = MU*(0.9*A). 

In the finite element grid each super-element is sub- 
divided N*N. The final elements are four-noded quadrilaterals 
using linear shape functions. For the transitional elements 
the additional nodes are found by using quadratic isoparametric 
shape functions for mapping.* The square grid on which the 
element mesh is based includes N*(2*N+1) nodes which are 
inside the structure. It is advantageous to retain these nodes 
in the numbering scheme. 

All of the nodal locations within the transition region 
can be derived from those in the sector from 0° to 45° by 
successive reflections. It follows that there are only 
(3*N+l)*N/2 distinct transition elements for which fluid "mass" 
and "stiffness" matrices need to be calculated. The remaining 
elements are identical squares and only a single additional 
pair of fluid matrices is required. 

A program (MGNEW) has been written to generate this 
element grid. Required input data are the structure radius A 

*Mapp ing with isoparametric shape functions is described in 
Ref . 3 . 

? 
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Fig. 1. Fluid region nomencl ature 
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Fig. 2. Super-element grid 
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and the mesh descriptors ML, MR, MU, and N. Output data in- 
clude nodal coordinates, numbers of nodes and elements, and 
global node numbers for each element. 

2 . 2 Fluid Program (DPLPOT) 

A program for tracking shock wave propagation, including 
cavitation, has been developed. Equations for the displacement 
potential formulation used are given in Ref. 1. 

2.21 Initial conditions 

A subroutine of the program determines initial values 
of the displacement potential and its first temporal derivative. 
These correspond to a uniform hydrostatic pressure plus a 
plane shock wave with step rise and exponential decay. The 
shock wave is travelling in the negative x direction with 
front at x = A at time t = 0. 

I 

j 

2.22 Boundary conditions 

The governing equation requires that values of the 
normal derivative of the displacement function be specified at 
boundary nodes. Subroutines are provided to impose a radiation 
boundary condition on the entry and exit faces. At the entry 
face explicit provision is made for the continuing entry of 
the incident shock wave. The analytic basis for these boundary 
conditions is given in Ref. 1. A radiation boundary condition 
is also imposed at the top face, again including explicit 
provision for the passage of the incident shock wave. In the 
absence of specific intervention the solution scheme enforces 
a zero value of the normal derivative at boundary points . 

I 

Accordingly, no boundary inputs are required for the plane 
of symmetry y = 0. 
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2.23 Temporary structure 

Initial testing of the fluid program DPLPOT has been 
done with a simplified boundary condition at the fluid- 
structure interface. This condition is that the boundary 
pressure remains equal to the initial hydrostatic pressure. 
Physically this corresponds to substituting an inflated 
cylindric bag for the structure. Such a "hole” in the fluid 
can readily be made to induce cavitation. 

2.24 Coefficient matrices 

Coefficient matrices for the fluid are calculated at 
the element level. The fluid "mass" matrix is lumped (diagonal) 
rather than consistent. Total storage requirements for element 
matrices are quite modest because of the limited number of 
different elements. The diagonal system mass matrix is 
assembled in a single vector whose length is equal to the 
number of nodes. No system "stiffness" matrix is assembled. 

2.25 Time integration 

Numerical integration in time is based on central 
difference formulas. The calculation is organized in such a 
way that the nonlinearity introduced by cavitation does not 
affect the fluid mass or stiffness matrices. 3ecause the 
integration algorithm is explicit and the mass matrix is 
diagonal, no matrix decomposition is required. As a result, 
the central processor time is directly proportional to the 
number of fluid nodes. Processing time is a linear function 
of the number of time steps. "Overhead” time requirements 
for mesh generation, coefficient matrix calculation and 
establishment of initial values are approximately equivalent 
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to 9 time steps. 

2 . 3 Fluid Program (VELPOT) 



The initial form of the program DPLPOT, described above, 
required a computationally expensive damping calculation. 

To avoid this, the alternative of using the velocity potential 
as the dependent variable was explored. Pertinent equations 
are given in Ref. 1. For this purpose it was necessary to 
revise four of the subroutines of DPLPOT. The alternate 
version, called VELPOT, reduced the central processor time by 
35%. Although VELPOT results closely resembled those from 
DPLPOT, they were not identical. Efforts to account for the 
discrepancies led to the development of two one-dimensional 
programs, DIS and VEE , which were made to yield identical 
results . 

Subsequent development of DPLPOT replaced the inefficien 
damping calculation. This improvement eliminated the computa- 
tional advantage of VELPOT and VELPOT was abandoned. 

2 . 4 Structural Program (STRUK2) 

A structural program has been under development by 
Jack T. Waller, a thesis student of mine. This program 
utilizes trigonometric series representations of the loading 
and displacements of an orthotropic shell in a state of plane 
strain. It is now being tested separately, but has not yet 
been coupled with the fluid. 

3. Results 

Results obtained thus far are relevant to determining 
when cavitation will occur, how big the cavity will become. 
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and how long it will persist. For this purpose, various 
combinations of shock pressure, hydrostatic pressure, and 
shock decay length have been tested. Effects of region size 
and shape, element size, length of time step, damping parameter, 
and total integration time have also been examined. 

3 . 1 Line Printer Plot 

The program provides for both tabular and graphical 
output, but the latter has been the more useful in this phase 
of the investigation. It is a line printer plot of fluid 
nodal pressures, repeated at selected time intervals. For 
the purpose of this plot the nodes of the structural and 
transition regions are mapped into the topologically equivalent 
rectangular array. 

Fig. 3 shows the arrangement of the line printer plot. 

The left-hand column is the plane of symmetry (x axis of Fig. 1) 
and the right-hand column is the "top" face. The entry face 
is the first row and the exit face is the last ro’v. Along the 
left-hand edge is a rectangle filled with X's, 4 columns wide 
and 7 rows high. The X's denote dummy nodes inside the 
structure. The rows immediately above and below, and the 
column to the right (all marked H) contain the 17 structural 
nodes . 

The pressure map of Fig. 3 is for hydrostatic pressure 
0.2 MPa and shock pressure 2.0 Mpa with a decay length of 15 m. 
The pressure ranges for the mapping characters are shown at 
the right. Since the map is for t = 0 there are no pressures 
below hydrostatic. 
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PRESSURE CODE 
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over 2.10 



Fig. 3. Initial pressure map 
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3.2 Cavitation Example 



Fig. 4 exhibits a sequence of pressure maps at 16 ms 
intervals based on the initial conditions of Fig. 3. Note 
that the cavity develops rapidly in the first two frames, 
shrinks in the next two, and is completely closed by 30 ms. 

The succeeding frames in which the cavity opens again represent 
the first instance in which this behavior was encountered. 

It has been established by a longer duration run that the 
cavity does not open again after 144 ms. 
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Fig. 4. Cavitation example 
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